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Abstract 

We report a direct lattice calculation of the K to tttt decay matrix elements for both the AI = 1/2 
and 3/2 amplitudes Aq and A^ on 2+1 flavor, domain wall fermion, 16 3 x 32 x 16 lattices. This 
is a complete calculation in which all contractions for the required ten, four-quark operators are 
evaluated, including the disconnected graphs in which no quark line connects the initial kaon and 
final two-pion states. These lattice operators are non-perturbatively renormalized using the Rome- 
Southampton method and the quadratic divergences are studied and removed. This is an important 
but notoriously difficult calculation, requiring high statistics on a large volume. In this paper we 
take a major step towards the computation of the physical K — > tttt amplitudes by performing a 
complete calculation at unphysical kinematics with pions of mass 422 MeV at rest in the kaon rest 
frame. With this simplification we are able to resolve Re(^4o) from zero for the first time, with a 
25% statistical error and can develop and evaluate methods for computing the complete, complex 
amplitude Aq, a calculation central to understanding the A = 1/2 rule and testing the standard 
model of CP violation in the kaon system. 

PACS numbers: 11.15.Ha, 12.38.Gc 14.40.Be 13.25.Es 
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I. INTRODUCTION 



The Cabibbo-Kobayashi-Maskawa (CKM) theory for the weak interactions of the quarks 
when combined with QCD provides a framework describing in complete detail all the prop- 
erties and interactions of the six quarks. This framework incorporates the most general 
assignment of masses and couplings and appears able to explain all observed phenomena 
in which these quarks participate. However, to date, the non-perturbative character of low 
energy QCD has obscured many of the consequences of the CKM theory. In particular, both 
the direct CP violation seen in K meson decay and the factor of 22.5 enhancement of the 
/ = 0, K — > Tin decay amplitude A relative to the 1 = 2 amplitude A 2 (the A J = 1/2 rule) 
lack a quantitative explanation. 

Wilson coefficients evaluated at a QCD scale of about 2 GeV represent the short distance 
physics and can be evaluated from the CKM theory using QCD and electro- weak perturba- 
tion theory. However, these factors explain only a factor of two enhancement of the 1 = 
amplitude 1 , [2 1 . The remaining enhancement must arise from the hadronic matrix elements 
which require non-perturbative treatment. 

Direct CP violation in kaon decays provides a critical test of the standard model's CKM 
mechanism of CP violation. While forty years of experimental effort have produced the 
measured result Re(e'/e) = 1.65(26) x 10 -3 with only a 16% error, there is no reliable 
theoretical calculation of this quantity based on the standard model. A previous lattice 
QCD calculation using 2+1 dynamical domain wall fermions failed to give a conclusive 
result because of the large systematic errors associated with the use of chiral perturbation 



theory at the scale o 
perturbation theory 



the kaon mass [4|. (However, there are on-going efforts using chiral 
5)].) Earlier quenched results jg, Q] are subject to this same difficulty 



together with uncontrolled uncertainties associated with quenching [8HlO]. 

A direct lattice calculation of K — > txtx decay is extremely important to provide an ex- 
planation for the AI =1/2 rule and to test the standard model of CP violation from first 
principles. This is an unusually difficult calculation because of the presence of disconnected 
graphs. However, with the continuing increase of available computing power and the devel- 
opment of improved algorithms, calculations with disconnected graphs are now no longer 
out of reach. In fact, our recent successful calculation of the masses and mixing of the 
rf and 77 mesons 11[ was carried out in part to develop and test the methods needed for 
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the calculation presented here. In this paper, we present a first direct calculation of the 
complete K° — > decay amplitude. At this stage, we work with the simplified kinematics 
of a threshold decay in which the kaon is at rest and decays into two pions each with zero 
momentum and with mass one-half that of the kaon. The calculation with this choice of 
kinematics still contains the main difficulties we need to overcome in order to be able to 
compute the physical K — > nir decay amplitudes; i.e. the presence of disconnected diagrams 
coupled with the need to subtract ultraviolet power divergences. However, as explained 
below, with the pions at rest we are able to generate sufficient statistics to explore how to 
handle these difficulties. We stress that at this simplified choice of kinematics, we compute 
the K — > 7T7T amplitudes directly and completely. 

In order to calculate the decay amplitudes, we perform a direct, brute force calculation 
of the required weak matrix elements. The isospin zero n — n final state implies the presence 
of disconnected graphs in correlation functions and makes the calculation very difficult. 
For these graphs, the noise does not decrease with increasing time separation between the 
source and sink, while the signal does. Therefore, substantial statistics are needed to get 
a clear signal. This difficulty is compounded by the presence of diagrams which diverge as 
1/a 2 as the continuum limit is approached (a is the lattice spacing). While these divergent 
amplitudes must vanish for a physical, on-shell decay they substantially degrade the signal 
to noise ratio even for an energy-conserving calculation such as this one. Studying the 
properties of the 1/a 2 terms and learning how to successfully subtract them is one of the 
important objectives of this calculation. The chiral symmetry needed to control operator 
mixing is provided by our use of domain wall fermions. 

Recognizing the difficulty of this problem, we choose to perform this first calculation on 
a lattice which is relatively small compared to those used in other recent work and to use a 
somewhat heavy pion mass (m^ 421 MeV) so we can more easily collect large statistics. 
We concentrate on exploring and reducing the statistical uncertainty since the primary goal 
of this work is to extract a clear signal for these amplitudes. Therefore, the quoted errors 
on our results are statistical only. 

The main objective of this paper is to calculate the AI = 1/2 decay amplitude A . A 
calculation of the AI = 3/2 part is included here for comparison and completeness. A much 
more physical calculation of this A J = 3/2 amplitude alone can be found in 12]. In the case 
of the 1 = 2 final state no disconnected diagrams appear, there are no divergent eye diagrams 
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and isospin conservation requires that four valence quark propagators must join the kaon 
and weak operator with the operators creating the two final-state pions. This allows physical 
kinematics with non-zero final momenta to be achieved by imposing anti-periodic boundary 
conditions on one species of valence quark jl^, [ijj]. As a result, the preliminary calculation 



of A 2 reported in Ref. [12| is performed at almost physical kinematics on a lattice of spatial 
size 4.5 fm and determines complex A 2 with controlled errors of 0(10%). The present work is 
intended as the first step toward an equally physical but much more challenging calculation 
of A Q . 

While we do not employ physical kinematics, the final results for the complex ampli- 
tudes A and A 2 presented in this paper are otherwise physical. In particular, we use 



Rome-Southampton methods 



15] to change the normalization of our bare lattice four-quark 



operators to that of the RI/MOM scheme. A second conversion to the MS scheme is then 



performed using the recent results of Ref. 16) . Finally these MS-normalized matrix elements 



3, 



determined in this same scheme, 



are combined with the appropriate Wilson coefficients 
to obtain our results for Aq and A 2 . Because of our unphysical, threshold kinematics and fo- 
cus on controlling the statistical errors associated with the disconnected diagrams, we do not 
estimate the size of possible systematic errors. Similarly we do not include the systematic 
or statistical errors associated with the Rome- Southampton renormalization factors, both 
of which could be made substantially than our statistical errors when required. 

This paper is organized as follows. We first summarize our computational setup, including 
our strategy to collect large statistics. Next we discuss our results for -n — -n scattering which 
are a by-product of the necessary characterization of the operator creating the tt — n final 
state and are also needed to evaluate the Lellouch-Luscher, finite- volume correction 18 ]. 
After a section giving the details of the K° — > tttt contractions, we provide our numerical 
results for the K° — > tttt decay amplitudes for both the AI = 3/2 and 1/2 channels. The 
details of the operator renormalization required by the Wilson coefficients which we use are 
presented in Appendix A. Finally we present our conclusions and discuss future prospects. 

II. COMPUTATIONAL DETAILS 

Our calculation uses the Iwasaki gauge action with /3 = 2.13 and 2+1 flavors of domain 
wall fermions (DWF). While the computational costs of DWF are much greater than those 
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of Wilson or staggered fermions, as has been shown in earlier papers [6|, |7J, L3J1 |20| , accurate 
chiral symmetry at short distances is critical to avoid extensive operator mixing, which 
would make the lattice treatment of AS* = 1 processes much more difficult. 

We use a single lattice ensemble with space-time volume 16 3 x 32, a fifth-dimensional 
extent of L s = 16 and light and strange quark masses of mi = 0.01, m s = 0.032, respec- 
tively. This ensemble is similar to the mi = 0.01 ensemble reported in Ref. 21] except we 
use the improved RHMC-II algorithm of Ref. [22] and a more physical value for the strange 
quark mass. The inverse lattice spacing for these input parameters was determined to be 
1.73(3)GeV and the residual mass is m res = 0.00308(4) [22j. The total number of configura- 
tions we used is 800, each separated by 10 time units. We initially generated an ensemble 
one-half of this size. When our analysis showed a non-zero result for Rev4 , we then doubled 
the size of the ensemble to assure ourselves that the result was trustworthy and to reduce the 
resulting error. We have performed the analysis described below both by treating the results 
from each configuration as independent and by grouping them into blocks. The resulting 
statistical errors are independent of block size suggesting that the individual configurations 
are essentially uncorrelated for our observables. 

We use anti-periodic boundary conditions in the time direction, and periodic boundary 
conditions in the space directions for the Dirac operator. The propagators (inverses of the 
Dirac operator) are calculated using a Coulomb gauge fixed wall source (used for meson 
propagators) and a random wall source (used to calculate the loops in the type3 and typei 
graphs shown in Figs. [5] and [6] below) for each of the 32 time slices in our lattice volume. For 
each time slice and source type, twelve inversions are required corresponding to the possible 
3 color and 4 spin choices for the source. Thus, all together we carry out 768 inversions 
for each quark mass on a given configuration. As will be shown below, this large number 
of inversions, performed on 800 configurations, provides the substantial statistics needed to 
resolve the real part of the 1 = amplitude A with 25% accuracy. 

The situation described above in which 768 Dirac propagators must be computed on a 
single gauge background is an excellent candidate for the use of deflation techniques. The 
overhead associated with determining a set of low eigenmodes of this single Dirac operator 
can be effectively amortized over the many inversions in which those low modes can be used. 
Our mi = 0.01, light quark inversions are accelerated by a factor of 2-3 by using exact, 
low-mode deflation 23j in which we compute the Dirac eigenvectors with the smallest 35 
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TABLE I: Masses of pion and kaons and energies of the two-pion states. Here the subscript 1 = 
or 2 on the ir — tt energy, EJ 77 , labels the isospin of the state and B™' represents the isospin zero, 
two-pion energy obtained when the disconnected graph V is ignored. The superscript (0), (1) or 
(2) on the kaon mass distinguishes our three choices of valence strange quark mass, m s = 0.066, 
0.099 and 0.165 respectively. 



m w E™ 




EY m§ m£ } 


(2) 

rrv K ' 


0.24373(47) 0.443(13) 


0.4393(41) 


0.5066(11) 0.42599(42) 0.50729(44) 


0.64540(49) 



eigenvalues and limit the conjugate gradient inversion to the remaining orthogonal subspace. 

In order to obtain energy- conserving K° — > tttt decay amplitudes, the mass of the valence 
strange quark in the kaon is assigned a value different from that appearing in the fermion 
determinant used to generate the ensembles, i.e. the strange quark is partially quenched. 
Since the mass of the dynamical strange quark is expected to have a small effect on ampli- 



tudes of the sort considered here 
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24J , this use of partial quenching is appropriate for the 
purposes of this paper. Valence strange quark masses are chosen to be m s = 0.066, 0.099 
and 0.165, which are labeled 0, 1 and 2 respectively. The resulting kaon masses are shown 
in Tab. HI In the following section we will see that by using these values for m s we can 
interpolate to energy- conserving decay kinematics for both the 1 = 2 and 1 = channels. 

III. TWO-PION SCATTERING 

The 7r — 7i scattering calculation requires 4 contractions which we have labeled direct 



(D), cross (C), rectangle (R), and vacuum (V) as in Ref. [25| and which are shown in Fig. [TJ 
For convenience, the minus sign arising from the number of fermion loops is not included 
in the definition of these contractions. The vacuum contraction should be accompanied by 
a vacuum subtraction. These contractions can be calculated in terms of the light quark 
propagator L(t snk , t STC ) for a Coulomb gauge fixed wall source located at the time t src and a 
similar wall sink located at t sn ^. The resulting complete vacuum amplitude, including the 
vacuum subtraction, is given by 

1 31 f 

V ® = (tr[L(t'X)L(t',tyML(t + t',t + t')L(t + t',t + ty]) (1) 

t'=0 I 
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FIG. 1: The four diagrams which contribute to tt — ir scattering: direct (D), cross (C), rectangle 
(R), and vacuum (V), arranged from the left top to right bottom. 

- (tr[L(t', t')L(t', t') f ] ) (tr [Lit + t',t + t')L(t + t',t + t') 1 "] ) J , 

where the indicated traces are taken over spin and color, the hermiticity properties of the 
domain wall propagator have been used to eliminate factors of 7 5 and we are explicitly 
combining the results from each of the 32 time slices. 

Our results for each of these four types of contractions are shown in the left panel of 
Fig. Notice that the disconnected (vacuum) graph has an almost constant error with 
increasing time separation between the source and sink, so it appears to have an increasing 
error bar in the log plot, while the signal decreases exponentially. 

These four types of correlators can be combined to construct physical correlation functions 
for two-pion states with definite isospin: 

(or(t+tyor(t')) = 2(d(*)-c(*)) (2) 

(0™(t + tyO™(t)) = 2D(t) + C(t) - 6R(t) + 3V(t). (3) 

Here the operator Oj n (t) creates a two-pion state with total isospin / and ^-component of 
isospin I z = using two quark and two anti-quark wall-sources located at the time-slice t. 
As in Eq. [1] we will average over all 32 possible values of common time displacement t' to 
improve statistics. 

The two-pion correlation functions for isospin / and I z = are fit with a functional form 
Corrr(t) = iV|{exp(— EJ^t) + exp(— Ej n (T — £)) + Cj}, where the constant Cj comes from 
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FIG. 2: Left: Results for the four types of contractions, direct (D), cross (C), rectangle (R), 
and vacuum(V) represented by the graphs in Fig. [TJ Right: Effective mass plots for correlation 
functions for states with isospin two (I2), isospin zero (Iq), isospin zero without the disconnected 
graph (7q) and twice the pion effective mass (2m n ). 

the case in which the two pions propagate in opposite time directions. The fitted energies 
are summarized in Tab. [H In order to see clearly the effect of the disconnected graph, we 
also perform the calculation for the 1 = channel without the disconnected graphs. This 
result is given in Tab. [I] with a label with an additional prime (/) symbol. The resulting 
effective mass plots for each case are shown in the right panel of Fig. |2j For comparison, a 
plot of twice the pion effective mass is also shown. This figure clearly demonstrates that the 
two-pion interaction is attractive in the 1 = channel with the finite volume, I = tt — ix 
energy Eq* lower than 2m n . In contrast, the 1 = 2 channel is repulsive with larger than 
2m^. The fitted parameters iV} r7r and E] n will be used to extract weak matrix elements from 
the K — > mi correlation functions discussed below in which these same operators 0* n (t) 
are used to construct the two-pion states. 

IV. CONTRACTIONS FOR K° -> vrvr DECAYS 



The effective weak Hamiltonian describing K° — > txtx decay including the u, d, and s 
flavors as dynamical variables is 

10 
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FIG. 3: Diagrams representing the eight K° — > itit contractions of typel, where Ty-^ = 7^(1 ±75). 
The black dot indicates a 75 matrix, which is present in each operator creating or destroying a 
pseudoscalar meson. 

Throughout this paper we follow the conventions and notation of Ref. [6|. In Eq. H]the Qi 
are the ten conventional four-quark operators, and tji are the Wilson coefficients, and r 
represents a combination of CKM matrix elements: r = —V t * s V t d/V u dV* s . To calculate the 
decay amplitudes A2 and Aq, we need to calculate the matrix elements (Tnr\Qi\K°) on the 
lattice. 

We list all of the possible contractions contributing to the matrix elements (7i7i\Qi\K°) 
in Figs. |3j|6j There are 48 different contractions which are labeled by circled numbers 
ranging from 1 to 48, and grouped into four categories labeled as typel, type2, type3, and 
typei according to their topology. Once we have calculated all of these contractions, the 
correlation functions {Oj n (tT r )Qi(t op )K (tK)) are then obtained as combinations of these 
contractions. In order to simplify the following formulae, we use the amplitude Aj ti (t n , t, tx) 
to represent three point function [O] 71 (t n )Qi(t op )K(t K )) . Using this notation, the 1 = 2 
amplitudes can be written, 



A2,l(t 7T , t op , tjc) 



gl©"©} 



(5a) 
(5b) 




FIG. 4: Diagrams for the eight type2 K° — > -kit contractions. 

A 2 ,3(tn,t op ,t K ) = (5c) 

A 2 ,a(U, t op , t K ) = (5d) 

A 2 ,5{U,t p,t K ) = (5e) 

A 2 fi(tn, top, ti<) = (5f) 

A 2 j(t 7T ,t op ,t K ) = i^{®-®} (5g) 

A 2 ${t n ,t ov ,t K ) = i\J^{®-®} (5h) 

A 2 ${t n ,t ov ,t K ) = «y|{0-©} (5i) 

^2,io(^,to P ,tx) = ^y|{©-©} (5j) 

and in the 1=0 case, 

A ,i{t^to P ,t K ) = *-J={-©-2-© + 3-® + 3-@-3-©} (6a) 

A , 2 (^,top,^) = *-J={-©-2-© + 3-@ + 3-@-3-©} (6b) 

A) )3 (^, io P , f *) = iVS{-® + 2- ®- @ + 2- @ + © (6c) 

— © — © — 2-© — © + © + ©} 

10 
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where the factor i comes from our definition of the interpolation operator for the mesons, 
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e.g. K° = i(dj 5 s). 

A few notes about the contractions shown in the Figs. [3] -[6] may be useful: 

1. The contractions identified by circled numbers do not carry the minus sign required 
when there is an odd number of fermion loops. Instead, the signs are included explicitly 
in Eqs. [5] and El 

2. The routing of the solid line indicates spin contraction while that of the dashed line 
indicates the contraction of color indices. If there is no dashed line, then solid line 
indicates connections implied by the trace over both color and spin indices. (This will 
be explained in more detail below.) 

3. A line represents a light quark propagator if it is not explicitly labeled with V. Up 
and down quarks and particular flavors of pion are not distinguished in Figs. [3] - [61 
Instead these specific contractions of strange and light quark propagators are combined 
in Eqs. [5] and [6] to give the 1 = 2 and 1 = amplitudes directly. 

4. Using Fierz symmetry, it can be shown that there are 12 identities among these con- 
tractions: 

© = -0, © = -©, © = -©, © = -©, (7a) 
@ = -@, © = -@, © = -©, © = -©, (7b) 
© = -©, © = -©, © = -©, © = -©■ (7c) 

A consequence of these identities is that Eq. E] is consistent with only seven of the ten 
operators Qi being linearly independent and with the three usual relations: 

Qw -Qg = Qi-Qs (8a) 
Q4-Q3 = Q2- Qi (8b) 
2Q 9 = 3Qj - Q 3 . (8c) 

5. Based on charge conjugation symmetry and 7 s hermiticity, the gauge field average of 
each of these contractions is real. 
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6. The loop contractions of type3 and typeA are calculated using the Gaussian, stochastic 
wall sources described in Sec. [TO 

In order to make our approach more explicit, we will discuss some examples. First 
consider the two contractions of typel identified as © and © and shown in the top half of 



where tx is the time of the kaon wall source, t w the time at which the two pions are absorbed 
and x op = (x op ,t op ) the location of the weak operator. The function L(x sink , i src ) is the 
light quark propagator, a 12 x 12 spin-color matrix, while S(x B i^k,t SIC ) is the strange quark 
propagator. The hermitian conjugation operation, f, operates on these 12 x 12 matrices. We 
use Tr c to indicate a color trace, Tr s a spin trace, and Tr, with no subscript, stands for both a 
spin and color trace. We have also used the 7 5 hermiticity of the quark propagators to realize 
the combination of quark propagators given in Eqs. [9] and [TOj, allowing both contractions 
to be constructed from light and strange propagators computed using Coulomb gauge fixed 
wall sources located only at the times t n and tx- Note the sum over the spatial components 
of the sink x n creates a symmetrical wall sink provided that the appropriate Coulomb gauge 
transformation matrix has been applied to the sink color index of this propagator to duplicate 
the Coulomb gauge transformation that was used to create the Coulomb gauge fixed wall 
source. We will sum over the spatial location, x op , of the weak operator, to project onto 
zero spatial momentum and improve statistics. Below we will show results as a function of 
the separations between t n , t op and tx- 

As a third example, which illustrates the use of random wall sources, consider contraction 
@ shown in Fig. [51 Using the notation introduced above, this contraction is given by 



Fig. EJ 




(10) 



(9) 



Tr s |7 M (l - 7 5 )L(x op ,t 7r )75 ^ L((x^, t n ), t K ) S(x op ,tx) j } 
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© = Tr{ 7M (l+75)^(Wo P )}^o P r (11) 



Here 77 (x) is the value of the complex, Gaussian random wall source at the space-time 
position x, while L H (x sink , t SIC ) is the propagator whose source is 7](x)5(xq — t STC ). The Dirac 
delta function 5(xo — t SIC ) restricts the source to the time plane t = t src . In the usual way, 
the average over the random source rj(x) which accompanies the configuration average, will 



op; 



set to zero all terms in which the source and sink positions for the propagator L (x Q p, t 
in Eq. [11] differ, giving us the contraction implied by the closed loop in the top left panel of 
Fig. [5j By using 32 separate propagators each with a random source non-zero on only one 
of our 32 time slices we obtain more statistically accurate results than would result from a 
single random source spread over all times. 

An important objective of this calculation is to learn how to accurately evaluate the 
quark loop integration that is present in type3 and type4 graphs and which contains a 1/a 2 , 
quadratically divergent component. As can be recognized from the structure of the diagrams, 
these divergent terms can be interpreted as arising from the mixing between the dimension- 
six operators Qi (for all i but 7 and 8) and a dimension-3 "mass" operator of the form 
S75G?. Such divergent terms are expected and do not represent a breakdown of the standard 
effective Hamiltonian written in Eq. HJ In fact, given the good chiral symmetry of domain 
wall fermions all other operators with dimension less than six which might potentially mix 
with those in Eq. H] will vanish if the equations of motion are imposed. Therefore these 
operators cannot contribute to the Green's functions evaluated in Eqs. |5] and [6] where the 
operators in Hw are separated in space-time from those operators creating the K meson and 
destroying the ir mesons, a circumstance in which the equations of motion can be applied. 

The problematic operator s'y^d is not explictly removed from the effective Hamiltonian 
because, again using the equations of motion, sj^d can be written as the divergence of an 
axial current and hence will vanish in the physical case where the weak operator Hw carries 
no four-momentum and is evaluated between on-shell states. While we can explicitly sum 
the effective Hamiltonian density H-w over space to ensure Hw carries no spatial momentum, 
to ensure that no energy is transferred we must arrange that the kaon mass and two-pion 
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energy are equal. We may achieve this condition, at least approximately, but there will be 
contributions from heavier states, which are normally exponentially suppressed, but which 
will violate energy conservation and hence will be enhanced by this divergent sj^d term. 

Since S'j^d will not contribute to the physical, energy-conserving K — > tttt amplitude, 
there is no theoretical requirement that it be removed. The coefficient of this sj^d piece 
is both regulator dependent and irrelevant. The contribution of these terms in a lattice 
calculation of K — > tttt decay amplitudes will ultimately vanish as the equality of the initial 
and final energies is made more precise and as increased time separations are achieved. 
However, the unphysical effects of this S75G? mixing are much more easily suppressed by 
reducing the size of this irrelevant term than by dramatically increasing the lattice size and 
collecting the substantially increased statistics required to work at large time separations. 

A direct way to remove this 1/a 2 enhancement is to explicitly subtract an aiSj^d term 
from each of the relevant operators Qi where the coefficient a, can be fixed by imposing the 
condition: 

{0\Qt - aiS l5 d\K) = 0, (12) 

a condition that is typically required in the chiral perturbation theory for K — > tttt jf|. Of 
course, this arbitrary condition will leave a finite, regulator-dependent s'j^d piece behind in 
the subtracted operator Qi — afsj^d. However, this unphysical piece will not contribute to 
the energy-conserving amplitude being evaluated. Since it is no longer l/a 2 -enhanced its 
effects on our calculation will be similar to those of the many other energy non-conserving 
terms which we must suppress by choosing equal energy K and tttt states and using sufficient 
large time separation to suppress the contributions of excited states. 
Following Eq. [12] we will choose the coefficient a, from the ratio 

(Note, with this definition the coefficient ai is proportional to the difference of the strange 
and light quark masses.) Thus, we will improve the accuracy when calculating graphs of 
type?) and type4 by including an explicit subtraction term for those operators Qi where 
mixing with S75G? is permitted by the symmetries (all but Q-j and Q%): 

(0™(t n )Q t (t op )K°(t K )) sub = (OQ 7T (t 7r )Qi(t op )K°(t K )) - a t <O ^(^)s 7 5rf(top)^ (^)> • 

(14) 
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mix?> mixA 



FIG. 7: Diagrams showing the contractions needed to evaluate the subtraction terms. These are 
labeled mix3 and mixA and constructed from the type3 and typeA contractions by replacing the 
operator Qi and fermion loop with the vertex S750L 

We should recognize that there is a second, divergent, parity-even operator sd which mixes 
with our operators Qi. However, we choose to neglect this effect because parity symmetry 
prevents it from contributing to either the K — > tttt or K — > |0) correlation functions being 
evaluated here. 

The amplitude (OQ 7T (t 7T )~s'-f 5 d(t op )K (tK)) includes two contractions, one connected and 
one disconnected as shown in Fig. [TJ These terms, which arise from the mixing of the 
operators Qi with sj^d, are labeled mix3 and mixA. To better visualize the contributions 
from different types of contractions, we can write the right hand side of Eq. [TH symbolically 

as 

typel + type2 + type3 + typeA — a ■ (mix3 + mixA) 
= typel + typel + sub3 + subA, (15) 

where sub3 = type3 — a ■ mix3 and subA = typeA — a ■ mixA. Note, here and in later 
discussions we refer to the term being subtracted as "mix" and the final difference as the 
subtracted amplitude "sub". 



V. K° 7T7T AI = 3/2 AMPLITUDE 

As Eqs. [5] and [7a| show, the AI = 3/2 K° — > 2n decay amplitude includes only typel 
contractions and four of the correlation functions are related 

3 3 

A 2 ,io = A 2t9 = -A 2A = -A %2 . (16) 



17 



1 .8e+08 



1.6e+08 - 



1.4e+08 - 



Q1 1 — 0— 
Q7/12 ^ 
Q8/48 



& A -& — a_ 



1.2e+08 - 



1e+08 



8e+07 



* Q o a> tp — 



3e+07 



2.5e+07 



2e+07 - 



Q1 >— e- 
Q7/12 



Q8/48 



A 



o ® 9 



A A © 



1 .5e+07 



1e+07 



6 8 10 12 



2 4 



10 12 14 16 



A = 12 



A = 16 



FIG. 8: Plots of the AI = 3/2 K° — > it — tt correlation functions for kaon source and tt — tt sink 
separations of A = 12 (left panel) and 16 (right panel). The x-axis gives the time t specifying the 
time slice over which the operator, Qi(x,t), i = 1, 7, 8, is averaged. The results for the operator 
Qi are divided by 12, and those for Q$ by 48 to allow the results to be shown in the same graph. 
The correlators C2,j(A,i) are fit using the A = 12 data with a fitting range 5 < t < 7. The 
resulting constants are shown as horizontal lines in both the A = 12 and 16 graphs. We can see 
that the A = 16 data are consistent with those from A = 12, but receive large contributions from 
the around-the- world paths. 

Therefore, we need only to calculate ^2,1, ^7 and A2 t &. The corresponding three correlation 
functions, C2,i(A,t) for % = 1, 7 and 8, with the choice of for the kaon mass, are shown 
in Fig. [SI Here we exploit our propagator calculation for sources on each of the 32 time 
slices to compute C^A, t) from an average over all 32 source positions: 

1 31 

C 2 ,i(A, t) = —J2 A ^ =t' + A,t op = t + t', t K = t'). (17) 



t'=0 



In Fig. E] we plot C 2|i (A,t) for < t < A at fixed A = 12 or 16. Table [D shows that m$ is 
almost equal to the energy of / = 2, tt — tt state, so the 3-point correlation function C2.i(A, t) 
should be approximately independent of t in the central region where the time coordinate 
of the operator is far from both the kaon and the two-pion sources, <C t <C A. 
We fit the correlators C2,i(A,£) using a single free parameter Mf^ 2 ' l£lt : 



C 2 ,i(A,t) = Mf^N^NKe-^e"^*-^, 



(18) 
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TABLE II: Results for the lattice AI = 3/2, K — > tttt transition amplitudes obtained from fitting 
the 3-point correlation functions to the functional form given in Eq. [18] for the six operators with 
AI = 3/2 components. The second column gives the lattice matrix elements M 3 ^ 2 ' ldt (xl0~ 2 ) while 
the third and fourth column give their contributions to the real and imaginary parts of 



M 3/2,iat (xl0 _ 2) Re(^ 2 )(GeV) Im(A 2 )(GeV) 



1 


0.4892(16) 


-1.737(ll)e-08 





2 


= Mi 


6.665(42)e-08 





7 


6.080(18) 


2.422(16)e-ll 


4.070(26)e-14 


8 


21.26(6) 


-1.979(13)e-10 


-9.646(61)e-12 


9 


=1.5Mi 


-7.917(50)e-15 


5.185(24)e-13 


10 


=1.5Mi 


6.103(38)e-12 


-1.448(9)e-13 


Total 




4.911(31)e-08 


-5.502(40)e-13 



where N K , mx and N n7T , E n7T are determined by fitting the kaon and two-pion correlators 
respectively: 

31 

— ( K (t + 0^(0) = N K {e~ mKt + e- m ^ T ~^) (19) 
62 f =o 

1 31 

— V] (OF(t + t')0™{t')) = (e- E ^ + e - E ™ (T - l) + C) . (20) 
32 

t'=0 

The constant C arises when the two pions join the source at t' and sink ait + t' by traveling 
in opposite time directions as discussed below. The fitted results for the matrix elements 

M 3/2,lat from A = U 

are listed in Tab. [II] in lattice units. 
Figure El shows that for the operators Q-j and Q$ the larger separation, A = 16, between 
the kaon source and tt—tt sink gives a much shorter plateau region than the case A = 12. This 
behavior is inconsistent with the usual expectation that it is the contributions from excited 
states of the kaon and pion, contributions which should be suppressed for larger A, that cause 
the poor plateau. An alternative, consistent explanation attributes the shortened plateau 
region seen for A = 16 to the 'around-the- world' effect. This is the contribution to the 
correlation function in which the two-pion interpolating operator at the sink annihilates one 
pion and creates another (instead of annihilating two pions as in the K — > nn contribution 
we are seeking) and the process at the weak operator is Kit — ► -n (instead of K — >■ nn). 
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FIG. 9: Diagrams showing the dominant around-the- world paths contributing to graphs of typel. 
The space-time region between the kaon wall source at tx and its periodic recurrence at tx + T is 
shown, where T = 32 is the extent of the periodic lattice in the time direction. For this around- 
the-world path, one pion travels directly from the pion wall source at to the weak operator, 
represented by the grey dot at t op . However, the second pion propagates in the other direction in 
time, passes through the periodic boundary and combines with the kaon before reaching the weak 
operator at t op . 

While one pion travels from the weak operator to the tt — tt sink the second is created at the 
sink and travels forward in time, passing through the periodic boundary to reach the weak 
operator together with the kaon. The corresponding dominant path is shown in Fig. [9j The 
time dependence of this behavior can be estimated as 

~ Mf /2M N^N K e- m ^ T e' (E ^- m ^ (21) 

which is A independent but suppressed by the factor exp(— m n T), where N n is the analogue 
of Nk for the case of single pion production and T = 32 is the temporal extent of the lattice. 
In contrast, the physical contribution in Eq. [18] is suppressed by exp(— E^A). Thus, the 
second, standard term falls with increasing A and the two factors are of similar size when 
A = T/2. Therefore, we should expect to see a large contamination from such around- 
the-world effects in the A = 16 case, consistent with Fig. [HJ In both panels of that figure, 
we plot as three horizontal lines the fitted result from A = 12 for the three amplitudes 
^■3/2,iat ]\f^]\f K ex p _ Ai^ for i — 1, 7 and 8. The agreement between these lines and the 
short plateaus seen in the right-hand, A = 16 panel indicates consistency between these two 
values of A. 
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Additional evidence supporting this explanation for the short plateau in the case of 
A = 16 can be obtained by examining the explicit dependence on t given by Eq. [21] for 
the around-the-world contribution. Examining the exponential decay with t in the A = 16 
correlators plotted in the right panel of Fig. [HI for operators Q7 and Q$ we find a value for 
Ekk — iti-it varying between 0.4 and 0.5 depending on the choice of fit range. A more accurate 
value of 0.498(2) can be obtained by fitting the corresponding correlator for A = 20 and a 
fit range of 5 to 11. The strangeness-carrying state whose mass we have labeled Ekk can be 
formed from two quarks and must be parity even. Direct calculation of Ek-k from a scalar 
sd correlator yields Ek-k = 0.752(12) which is consistent with the sum of the result above, 
Ekk — m-n = 0.498(2), and the pion mass m n = 0.2437(5). (This energy difference is also 
close to the kaon mass m$ = 0.50729 given in Tab. HI) Thus, the time dependence expected 
from the around-the-world path is quite consistent with that seen in Fig. [HJ 

We conclude that it is important to increase the lattice extent in the time direction 
both to suppress this around-the-world effect and to permit the use of a larger source-sink 
separation giving a longer plateau. We will return to discussion of the around-the-world 
effect below for the A J = 1/2 kaon decay where it creates even greater difficulties. However, 
here we can begin to appreciate the severity of this effect in the K° — > system for our 
temporal lattice extent of 32, given our values of the lattice spacing and meson masses. 

The Wilson coefficients and operators which appear in Eq. H] are typically expressed in 
the MS scheme. Thus, we must change the normalization of our lattice operators Qi to that 
of the MS scheme. We begin by converting our bare lattice operators into the regularization 
invariant momentum (RI/MOM) scheme of Ref. [jjj]. Here we use the earlier results of 
Ref. (26 ] which were obtained for the present lattice action using the methods of Ref. jg]. 
In this previous work off-shell, Landau-gauge-fixed Green's functions containing the lattice 
operators Qi are evaluated at specific external momenta characterized by an energy scale /z. 
These results determine a renormalization matrix Z RI (/i, a) which can be used to convert 
the lattice normalization into that of the RI scheme: 

Q R V)* = E^r R W)Q;-- (22) 

As explained in Appendix A, these equalities hold only when the operators appear in phys- 
ical matrix elements. The indices i and j take on seven values corresponding to the seven 
independent operators in what will be called the chiral basis. (The primes in this equation 
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indicate lattice operators denned in that basis.) This is referred to as nonperturbative renor- 
malization (NPR) because the matrix Z]| t ~ ) ' RI (/i, a) is computed using a lattice evaluation 
of off-shell Green's functions and perturbation theory is not used. 

Next these Q R1 (fi)i operators are converted to the MS scheme in which the Wilson coeffi- 
cients are evaluated by applying a conversion matrix i? RI_s>MS discussed in detail in Ref. [if]]. 
Finally the matrix elements of these MS operators are combined with the Wilson coefficients 
obtained in the MS scheme 1?| using the scale p — 2.15 GeV to determine the results given 
later in this section for the A J = 3/2 amplitude A 2 and in the following section for the 
AI = 1/2 A . These procedures are described in greater detail in Appendix A. 

A good approximation to the infinite volume decay amplitude can be obtained by includ- 
ing the Lellouch-Luscher factor (F) 18] which relates the K — >■ tttt matrix element M of 
the effective weak Hamiltonian of Eq. H] calculated using finite volume states normalized to 
unity to the infinite volume amplitude A: \A\ 2 = F 2 M 2 where 



Here p is defined through = 2^Jm 2 + p 2 , q = Lp/2n and 5 2 (p) is the s-wave, 1 = 2, 
71 — 7r scattering phase shift for pion relative momentum p. The function <ft(q) is known 



analytically and given, for example, in Ref. 18]. The 1 = 2 phase shift ^(p) is determined 



from the measured two-pion energy E n7T = 0.443(13) given in Tab. [I] and the finite volume 



quantization condition 



27] 

0(?) + hip) = nn. (24) 



For our threshold case we set the integer n to zero and obtain ^(p) = —0.0849(43). Because 
of the small value of p we assume that ^(p) is a linear homogenous function of p and write 
hip) — pdhip) / dp, the quantity required in Eq. [23] and given in Tab. IHII (Equation |2"31 
differs by a factor of two from the expression given in the Lellouch-Luscher paper because of 
our different conventions for the decay amplitude A. With our conventions the experimental 
value of Re(A 2 ) = 1.48 x 10~ 8 GeV.) 

In the limit of non-interacting pions, the factor F becomes F f ^ ee = 2(2m 7r ) 2 mi<-Z/ 3 , which 
reflects the different normalization of states in a box and plane wave states in infinite volume. 
Results for F in this 1 = 2 case and the quantities used to determine it are given in Tab. IHII 
We should note that applying the finite volume correction of Eq. [23] gives us a finite- volume 
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TABLE III: The calculated quantities which appear in the Lellouch-Liischer factor F for 1 = 2. 
The corresponding factor for the case of non-interacting particles is F{ Tee = 31.42. The difference 
reflects the final two-pion scattering in a box. 



p 


9<t>(l) 
9 dq 


as{ P ) 
P dp 


F 


0.0690(13) 


0.221(10) 


-0.0849(43) 


26.01(18) 



corrected amplitude for a AI = 3/2, K — > tttc decay that is slightly above threshold by the 
amount - 2m n = 33(1) MeV. 

We can now combine everything and calculate the K° — > -kit decay amplitudes, 



^2/0 



^ 10 7 

'^XEE [(^) + ry l ^))z^ m M^ M ] , (25) 



where the construction of the 10x7 renormalization matrix ^t^-MS j s explained in Appendix 
A. For later use we have written Eq. [25] in a way which is applicable for AI = 1/2 decays 
as well as for the AI = 3/2 transitions considered in this section. The results for the 
complex A J = 3/2 decay amplitude A 2 are summarized in Tab. IIV| including those for 
the other two, energy-non-conserving choices of kaon mass. Since rrS^ differs from the 
isospin-2 n — n energy by only 0.2 percent, we quote this case as our energy-conserving kaon 
decay amplitude. Therefore, in physical units, we obtain the energy- conserving AI = 3/2, 
K° 7T7T complex, threshold decay amplitude for mx = 877 MeV and m n = 422 MeV: 

Re(A 2 ) = 4.911(31) x 10~ 8 GeV (26) 
lm(A 2 ) = -0.5502(40) x 10~ 12 GeV. (27) 

This result for Re(^2) can be compared with the experimental value of 1.48 x 10~ 8 GeV 
given above. The larger result found in our calculation is likely explained by our unphysically 
heavy kaon and pions. 

VI. K° —7- 7T7T AI = 1/2 AMPLITUDE 

Following the prescription given by Eq. Owe have calculated all of the AI = 1/2 kaon 
decay correlation functions, 

1 31 

C ,i(A, t) = — 4),i&r = t' + A, t op = t + t K = t'), (28) 

t'=0 
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TABLE IV: The complex, K° — > tttt, AI = 3/2 decay amplitudes in units of GeV. 



m K 


Re(^ 2 )(xl0- 8 ) 


Im(^ 2 )(xl0- 12 ) 


(0) 

m K 

(i) 
m K 

(2) 

rrv K ' 


4.308(28) 
4.911(31) 
5.916(38) 


-0.5596(40) 
-0.5502(40) 
-0.5316(39) 



for each of the ten effective weak operators. In the calculation we treat each of these ten 
operators as independent and then verify that the identities shown in Eq. [8]are automatically 
satisfied. Figures [10] and [11] show two examples of the resulting correlation functions for the 
operators Q2 and Qq, in the case of the lightest kaon . Table [I] shows that the mass of 
this kaon is very close to the energy of the 1=0 two-pion state. Therefore, we expect to get 
a reasonably flat plateau when the operator is far from both the source and sink. 

Given this good agreement between the energies of the K and 71 — 71 states, we might 
expect that the unphysical, dimension three operator, s^d which mixes with the (8, 1) 
operators in Eq. H] and is itself a total divergence, will also give a negligible contribution to 
such an energy and momentum conserving matrix element. However, as can be seen from 
Figs. [TOT a) and fllTa). the matrix element of this term is large and the explicit subtraction 
described in Sec. [TV] is necessary. 

This difficulty is created by the combination of two phenomena. First the mixing coef- 
ficient which multiplies the s^d operator when it appears in our weak (8, 1) operators is 
large, of order (m s — mi) /a 2 . Second, in our lattice calculation the necessary energy con- 
serving kinematics (needed to insure that this total divergence does not contribute) is only 
approximately valid. The required equality of the spatial momenta of the kaon and 71 — 71 
states is assured by our summing the location of the weak vertex over a complete temporal 
hyperplane. On the other hand, the equality of the energies of the initial and final states 
results only if we have adjusted the kaon mass to approximately that of the two-pion state 
and chosen the time extents sufficiently large that other states with different energies have 
been suppressed. However, as can be seen in Figs. ITOT a) and [TTT a) the subtraction terms 
mix3 and mixA show strong dependence on the time at which they are evaluated. This 
implies that there are important contributions coming from initial and final states which 
have significantly different energies. One or both of these states is then not the intended K 
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(c) (d) 
FIG. 10: Plots showing the t dependence of the various contractions which contribute to the 
AI = 1/2 correlation function Co, 2 (A = 16, t) for the operator (a) Contractions of type3, the 
divergent mixing term mix3 that will be subtracted and the result after subtraction, sub3. (b) 
Contractions of type4, the divergent mixing term mixA that will be subtracted and the result after 
subtraction, subA. (c) Results for each of the four types of contraction after the needed subtractions 
have been performed, (d): Results for the complete Q2 correlation function Cb,2(A = 16, t) 
obtained by combining these four types of contractions. The solid points labeled Q% are the 
physical result while the open points labeled Q' 2 are obtained by omitting all the vacuum graphs, 
sub4. The solid and dotted horizontal lines indicate the corresponding fitting results and the time 
interval, 5 < t < 11 over which the fits are performed. 
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(c) (d) 
FIG. 11: The result for each type of contraction contributing to the 3-point correlation function 
Co, 6 (A = 16, t) for the operator Qq following the same conventions as in Fig. 1101 



or 7i — 7r state but instead an unwanted contribution which has been insufficiently suppressed 
by the time separations between source, weak operator and sink. 

Thus, instead of relying on large time extents and energy conserving kinematics to sup- 
press this unphysical, 0(l/a 2 ) term we must explicitly remove it. As explained in Sec. [IV 
this can be done by including an explicit subtraction which we fix by the requirement that 
the kaon to vacuum matrix element of the complete subtracted operator vanishes as in 
Eq. [12j Thus, we determine the divergent coefficient of this mixing term from the ratio 
®-i — {®\Qi\K) / (0|s7 5 <i|fQ and then perform the explicit subtraction of the resulting terms, 
labeled a, • mix3 and «j • mix4 in Figs. [10] and [TTJ 

Of course, the finite part of such a subtraction is not determined from first principles 
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and our choice, specified by Eq. [12] is arbitrary. Thus, we must rely on our identification 
of a plateau and the approximate energy conservation of our kinematics to make the arbi- 
trary part of this subtraction small, along with the other errors associated with evaluating 
the decay matrix element of interest between initial and final states with slightly different 
energies. 

We now examine the very visible time dependence in Figs. flOT a) and [TTT a) for both the 
original matrix elements and the subtraction terms in greater detail. As discussed above one 
might expect these divergent subtraction terms to contribute to excited state matrix elements 
in which the energies of the initial and final states are very different. Typical terms should 
be exponentially suppressed as the separation between the weak operator and the source or 
sink is increased, with the time behavior exp{ — (m* K — mx)t} or exp{ — (E* n — E n7T )(A — t)}, 
which ever is larger. (The * denotes an excited state.) However, by carefully examining 
the time behavior of the mix3 amplitude, we find that the time dependence, at least in the 
vicinity of the central region, is less rapid than might be expected from such excited states 
suggesting that it is probably not due primarily to contamination from excited states. 

We believe that the dominant, energy-nonconserving matrix elements which cause the 
significant time dependence in Figs. [10] and [11] arise from the around-the-world effects iden- 
tified and discussed in the previous AI = 3/2 section. In fact, for the reasons just discussed 
associated with divergent operator mixing, such around-the-world effects are a more serious 
problem in the AI = 1/2 case. The dominant around-the-world graphs are shown in Fig.[T2l 
An estimate of the time dependence of these graphs gives, 

< K°7r\Qi\7r > N^N K N v e- m * T e^^-™^ 

+ < 0\Qi\K°7T7T > N^N K N^e- mK ^ T -^ + ^-^ , (29) 

where the first term comes from the first two graphs of Fig. [12], while the second term comes 
from the third graph. (Recall that t = t op — t K and A = t w — t K ). Notice that these two 
terms involve amplitudes which are far from energy conserving and therefore contain large 
divergent contributions from mixing with the operator sj^d which will be removed only when 
combined with the corresponding around-the-world paths occuring in the mix3 contraction. 

We conclude that it is these around-the-world matrix elements which are the reason 
for the observed large divergent subtraction in the type3 graph. The largest divergent 
contribution is thus not the subtraction for the matrix element we are trying to evaluate, 
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FIG. 12: The dominant around-the-world paths contributing to graphs of type3. As in Fig. [9] 
we show the space-time region between the kaon source at t = tx and its periodic recurrence at 
t = tx + T. The gray circle represents the four quark operator Qi. For the first two graphs, one of 
the two pions created at the t = t n source travels directly to the operator Qi while the second pion 
travels in the other direction in time and reaches the kaon and weak operator by passing through 
the periodic lattice boundary. In the third diagram it is the kaon which travels in the opposite to 
the expected time direction. 

< 7T7r\Qi\K° >; rather, it is the divergent subtraction for the matrix elements < K°7r\Qi\7i > 
and < 0\Qi\K 7T7T > which arise from the around-the-world paths which are not sufficiently 
suppressed by our lattice size. Two important lessons can be learned from this analysis. 
First, it is important to perform an explicit subtraction of the divergent mixing with the 
operator s^^d. While this term will not contribute to the energy conserving matrix element 
of interest, in a Euclidean space lattice calculation there are in general, other, unwanted, 
energy non-conserving terms which may be uncomfortably large if this subtraction is not 
performed. Second it would be wise to work on a lattice with a much larger size T in time 
direction in order to suppress further the around-the-world terms which give such a large 
contribution in the present calculation. Using the average of propagators computed with 
periodic plus anti-periodic boundary conditions to effectively double the length in the time 
direction would be a good solution. 
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We should emphasize that these divergent, around-the- world contributions do not pose a 
fundamental difficulty. The largest part of these amplitudes are removed by the correspond- 
ing subtraction terms constructed from the operator s^d. The remaining finite contribu- 
tions from this and other around-the- world terms are suppressed by the factor exp(— m n T) 
or exp(— mk(T— A)). Fortunately, the large divergent subtraction also reduces the statistical 
errors substantially, especially for the typeA vacuum graphs, which indicates the expected 
strong correlation between the divergent part of the weak operator and the corresponding 
S75G? subtraction. Our results suggest that the separation of A = 16 gives a relatively longer 
plateau region, so we use that K — im time separation in the analysis below. 

The lattice matrix elements are determined by fitting the I = 1/2 correlators Cq(A,£) 
given in Eq. [28] using the fitting form: 

C 04 (A,t) = M^N^NKe-^e-^-^. (30) 

The fitted results for the weak, AI = 1/2 matrix elements of all ten operators are sum- 
marized in Tab. [V] To see the effects of the disconnected graph clearly, a second fit is 
performed to the amplitude from which the disconnected, type4 graphs have been omitted 
and the calculated results are shown with an additional / label, as in the earlier two-pion 
scattering section. 

The calculation of the AI = 1/2 decay amplitude Aq from the lattice matrix elements 
given in Tab. [V] is very similar to the AI = 3/2 case: the values of M± ' a are 
simply substituted in Eq. [25j However, the attractive character of the I = 0, 7r — it in- 
teraction and resulting negative value of p 2 makes the Lellouch-Liischer treatment of finite 
volume corrections inapplicable. For the repulsive 1 = 2 case, we could apply this treatment 
to obtain the decay amplitude for a two-pion final state which was slightly above threshold 
corresponding to the actual finite volume kinematics. In the present case there is no corre- 
sponding infinite- volume decay into two pions below threshold and an unphysical increase of 
to compensate for the finite volume 7r — rr attraction will introduce an 0(1/L 3 ) error in 
the decay amplitude of the same size as that which the Lellouch-Liischer treatment corrects. 
Thus, for this AI = 1/2 we do not include finite volume corrections and simply use the 
free-field value for the factor F in Eq. 1251 

While we believe that we cannot consistently apply the Lellouch-Liischer finite volume 
correction factor to improve our result for the I = 0, K — > tttt decay amplitude, we might 
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TABLE V: Fitted results for the weak, AI = 1/2 kaon decay matrix elements using the kaon 
mass rripr- The column M\ &t shows the complete result from each operator. The column M 4 - lat 



shows the result when the disconnected graphs are omitted while the 4th and 5th columns show 
the contributions of each operator the real and imaginary parts of the physical decay amplitude 
Aq. These results are obtained using a source-sink separation A = 16, and a fit range 5 < t < 11. 



i 


, .-l/2,lat/ 1n _2\ 

Mj (xlO J 


M i (xlO ) 


T"> (A \ ( \T\ 

Ke(74oj(CjeVJ 


Im(^o)(GeV) 


1 


-1.6(16) 


-1.10(37) 


7.6(64)e-08 





2 


1.52(61) 


1.92(15) 


2.86(97)e-07 





3 


-0.3(41) 


0.3(10) 


2.1(136)e-10 


l.l(76)e-12 


4 


2.7(33) 


3.32(78) 


4.2(44)e-09 


1.4(14)e-ll 


5 


-3.3(38) 


-6.81(86) 


3.1(53)e-10 


1.6(28)e-12 


6 


-7.8(48) 


-19.6(9) 


-5.6(33)e-09 


-3.3(20)e-ll 


7 


10.9(14) 


15.20(42) 


5.2(12)e-ll 


8.8(20)e-14 


8 


35.7(28) 


47.2(10) 


-3.66(28)e-10 


-1.79(14)e-12 


9 


-2.2(12) 


-1.79(29) 


3.1(15)e-14 


-2.01(96)e-12 


10 


0.9(12) 


1.24(29) 


1.2(ll)e-ll 


-2.7(27)e-13 


Total 






3.46(78)e-07 


-2.4(23)e-ll 



still be able to use the quantization condition of Eq. EH to determine the I = tt — n 
scattering phase shift Sq (p) . Even though Eq. |2U can be analytically continued to imaginary 
values of the momentum p, its application for large negative p 2 is uncertain since the function 
<p(q) becomes ill denned. In fact, our value of p 2 sits very close to a singular point of <f)(q). 
We believe this happens because the condition on the interaction range ^ < L/2 used to 
derive the quantization condition in Eq. [2H is not well satisfied for our small volume. This 
impediment to determining 8o(p) will naturally disappear once we work with lighter pions 
in a larger volume. 

The results for Re(A ) and Im(A ) are summarized in Tab. |VT]and the individual contri- 
bution from each of the operators is detailed in the last two columns of Tab. |Vl Within a 
large uncertainty Tab. IVl shows that the largest contribution to RefAi) comes from operator 
Q2, and that to Im(v4 ) from Q 6 as found, for example, in Refs. p, |7]. 

Since the choice for the kaon mass is not precisely equal to the energy of the I = 
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TABLE VI: Amplitudes for AI = 1/2 K° -> vrvr decay in units of GeV. The energy conserving 
amplitudes are obtained by a simple linear interpolation between =0.42599 and =0.50729 
to the energy of two-pion state. As in the previous tables, the / indicates results from which the 
disconnected graphs have been omitted. 

m K Re(A))(xlO- 8 ) Re(A^)(xHr 8 ) Im(A )(x lO" 12 ) Im(A^)(xlO- 12 ) 



m K {0) 


36.1(78) 


42.3(20) 


-21(21) 


-66.1(43) 


m K (l) 


45(10) 


48.8(24) 


-41(26) 


-74.6(47) 


m K {2) 


65(15) 


58.6(32) 


-69(39) 


-89.6(63) 


Energy conserving 


38.0(82) 


43.4(21) 


-25(22) 


-67.5(44) 



state, we carried out a simple linear interpolation between m K and m K to obtain an energy 
conserving matrix element, which is shown in the last row of Tab I VI I In terms of physical 
units, therefore, our full calculation gives the energy conserving, K° — > tttt, AI = 1/2, 
complex decay amplitude A for mx = 766 MeV and = 422 MeV: 

Re(A ) = 3.80(82) x l(T 7 GeV (31) 
Im(Ao) = -2.5(2.2) x lCr n GeV. (32) 

These complete results can be compared with those obtained when the disconnected graphs 
are neglected given in Tab. |VI] and the experimental value for Re(A)) = 3.3 x 10~ 7 GeV. 
As in the case of Ke(A2), our larger value is likely the result of our unphysically heavy kaon 
and pion. 

VII. DISCUSSION AND CONCLUSIONS 

Comparing the results of Re(A 2 ) in Tab. HVland Re(A ) in Tab. |VI| we find the AI = 1/2 
enhancement ratio Ke(A ) /R,e(A 2 ) to be roughly 7-9. This comparison is degraded by our 
threshold kinematics which, since the 1 = and 1 = 2 two-pion states have different energies 
in a finite volume, causes us to use a different kaon mass in the calculations of (^2) and 
(Ao) in order to have energy conserving decays in each case. These two energy conserving 
amplitudes have a ratio of 38.0/4.911 = 7.7, while if we ignore energy conservation and use 
the same value for kaon mass, the ratio becomes 45.0/4.911 = 9.2. Of course, both 
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estimates are far from the experimental ratio of 22.5 suggesting that our 422 MeV pion mass 
and small lattice volume are far from physical. 

For completeness, we also calculate the measure of direct CP violation, 



where u> = Re(A 2 ) /Re(A ) is the inverse of the AI = 1/2 enhancement factor. Using 
our kinematics, the kaon mass and substituting the experimental value for e, we get 
Re(e'/e) = (2.7 ± 2.6) x 10" 3 . If we instead use the experimental value for u, we get 
Re(e7e) = (1.11 ±0.91) x 1(T 3 . 

Our calculation is sufficiently far from physical kinematics, that it is not appropriate 
to compare these results with experiment. 1 Instead, our objective is to show how well our 
method performs. We have been able to calculate Re(A ), the key element needed to explain 
the AI = 1/2 rule, with a 25% statistical error. Comparing our results for Re(A ) obtained 
on sub-samples of N=100, 400 and all 800 configurations we find that the statistical errors 
on the quantities we measure do indeed scale as 1/ \^N. Therefore, we believe that our non- 
zero signal for Re(Ao) is real and that we could reduce this statistical error to 10 percent by 
quadrupling the size of our sample to 3200 configurations. It is interesting to note the results 
for primed (disconnected graphs omitted) and unprimed (all graphs included) quantities 
contributing to Re(A ) have similar values suggesting that the disconnected graphs, while 
contributing significantly to the statistical error, have an effect on the final result for Re(A ) 
at or below 25%. 

In contrast, the result for Im(Ao) has an 80% error. Thus, it is not clear whether the size 
of the result will survive a quadrupling of the sample with its statistical error reducing to 
a 40% error or whether the result itself will shrink, remaining statistically consistent with 
zero. Considering the substantial systematic errors associated with our small volume and 
the fact that our kinematics are far from the physical, we present this trial calculation as 
a guideline for future work and a proof of method rather than giving accurate numbers to 
compare with experiment. 

From our observation of the around-the- world effect, we conclude that it is important 

1 A further unphysical aspect of our kinematics is the inequality of the strange quark mass used in the 
fermion determinant and the self contractions appearing in the eye graphs (m s — 0.032) and strange 
quark masses used in the valence propagator of the K meson (m s = 0.066, 0.99 and 0.165). 




y/2\e\ [Re(A 2 ) Re(A ) 



u rim(A 2 ) Im(A)) 



(33) 
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to use the average of quark propagators obeying periodic and anti-periodic boundary con- 
ditions to extend the lattice size in the time direction. In addition, explicit subtraction of 
the divergent mixing term s^ 5 d is necessary even for kinematics which are literally energy 
conserving because the around-the-world path and possibly other excited state matrix ele- 
ments are far off shell and can be substantially enhanced by such a divergent contribution. 
Finally, future work should be done using a much larger lattice which can contain two pions 
without any worry about finite size effects. 

The focus of this paper is on developing techniques capable of yielding statistically mean- 
ingful results from the challenging lattice correlation functions involved in the amplitude Aq. 
However, there are other important problems that will also require careful attention if physi- 
cally meaningful results are to be obtained for this amplitude with an accuracy of better than 
20%. Two important issues are associated with operator mixing. As discussed in Appendix 
A, a proper treatment of the non-perturbative renormalization of the four independent (8, 1) 
four-quark operators requires that additional operators containing gluonic variables (some 
of which are not gauge invariant) be included. While including such operators is in principle 
possible and the subject of active research, controlling such mixing using RI/MOM methods 
offers significant challenges. 

A second problem is operator mixing induced by the residual chiral symmetry breaking of 
the D WF formulation. The mixing of such wrong- chirality operators should be suppressed 
by a factor of order m res . However, the K — )■ tctt matrix elements of the important (8, 1) four- 
quark operators are themselves suppressed by at least one power of m 2 K) a suppression that 
is absent from similar matrix elements of the induced, wrong-chirality operators. Therefore, 
such mixing has been ignored in this paper because its effect on the matrix elements of 
interest are expected to be of order m Tes /m s ~ 0.08, suggesting that these effects will be 
smaller than our 25% statistical errors. To perform a more accurate calculation in the 
future, these mixing effects may be further suppressed by adopting a gauge action with 
smaller residual chiral symmetry breaking. For example, this ratio reduces to 0.04 for the 



DSDR gauge action now being used in RBC/UKQCD simulations [28 



and to 0.023 for 



those ensembles with the smallest lattice spacing created to date using the Iwasaki gauge 
action 29] . When greater accuracy is required either an improved fermion action, larger L s 



or explicit subtraction of wrong-chirality mixing must be employed. 

As we move closer to the physical pion mass we must overcome a further important 
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difficulty: giving physical relative momentum to the two pions. This can be accomplished 
while keeping the two-pion state in which we are interested as the ground state, if the 
kaon is given non-zero spatial momentum relative to the lattice. In this case the lowest 
energy final state can be arranged to have one pion at rest while the other pion carries the 
kaon momentum, as in the AI = 3/2 calculation of Ref. jsoj. However, this requires the 
momentum carried by the initial kaon and final pion to be 739 MeV, which is 5.4 times 
larger than the physical pion mass. Such a large spatial momentum will likely make the 
calculation extremely noisy. For the AI = 3/2 calculation, it is possible to use anti-periodic 
boundary conditions in one or more spatial directions for one of the light quarks so that each 
pion necessarily carries the physical, 206 MeV momentum present in the actual decay while 



the kaon can be at rest 



13] . However, this approach cannot be used in the case of the 



1 = final state being studied here. Instead, the use of G-parity boundary conditions 
may be the solution to this problem. 
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Appendix A: Operator normalization 

In order to combine our lattice matrix elements with the Wilson coefficients describing 
the short- distance weak interaction physics responsible for K — > titt decay we must convert 
our lattice operators into those normalized according to that MS scheme in which the Wilson 
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coefficients are evaluated. We will discuss the details of this procedure in this appendix. 
The first step is converting the lattice operators into those normalized according to the 



RI/MOM scheme [15]. We follow the procedure of Ref. [6[ and make use of the fact that the 
ten operators which enter the conventional expression given in Eq. H] are linearly dependent 
and can be reduced to a set of seven independent operators, Q[, Q' 2 , Q' 3 , Q' 5 , Q' 6 , Q' 7 and 
Q' s defined in Eq. 172-175 Ref. jg]. These have been defined so that the resulting opera- 
tors belong to specific irreducible representations of SUl(3) x SUr(3). The operator Q[ 
transforms as a (27, 1). The four operators Q' 2 , Q' 3 , Q' 5 and Q' 6 all belong to the (8, 1) rep- 
resentation, while Q' r and Q' 8 each transform as an (8,8). Here (m,n) denotes the product 
of an m-dimensional irreducible representation of SUl(3) with an n-dimensional irreducible 
representation of SUr(3). We refer to the basis of these seven independent operators as 
the chiral basis. Because SUl(3) x SUr(3) is an exact symmetry of the large momentum, 
massless limit which our NPR calculation is intended to approximate, the mixing matrix 
Z lat_s>RI given in Eq. [22] which relates the lattice and Rl-normalized operators will be block 
diagonal, only connecting operators which belong to the same irreducible representation of 
SU L {3) x SU R {3). 

The RI/MOM conditions which define the operators Of 1 and determine the 7x7 matrix 
j£iat->M are imposed on the Green's functions: 2 

Gi(px,P2) f aM5 = II {/ (^i)af(x2)pQf(0)d 7 (x 3 )f s (x 4 )) e -^(^+^) e ^N+-4) 

(Al) 

evaluated for p\ = p\ = (pi — P2) 2 = H 2 - Here a, (3, 7 and 5 are spin and color indices. 
The fields d and / create a down quark and a quark of flavor / = u or d while s and / 
destroy a strange quark and a quark of flavor /. The RI/MOM conditions are imposed 
by removing the four external quark propagators from the amplitudes in Eq. IA11 and then 
contracting each of the resulting seven amputated Green's functions obtained from Eq 
with seven projectors {r^{, 5 }i<j<7. The matrix Z lat ^ Rl is then determined by requiring tl 
the resulting 49 quantities take their free field values, as is described in detail in Refs. 
and Q. 



2 While this equation agrees with Eqs. 143 and 152 of Ref. [6|, a different choice of momenta was actually 
used in that earlier reference. These two equations accurately describe the earlier kinematics only after 
one pair of the momenta p\ and P2 are exchanged: pi f> 
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TABLE VII: The renormalization matrix Z lat_>,RI /Zg in the seven operator chiral basis at the 
energy scale \i = 2.15 GeV. These values were obtained from Ref. [26] by performing an error 
weighted average of the values given in Tabs. 40, 41 and 42 (corresponding to bare quark masses 
of 0.01, 0.02 and 0.03) and inverting the resulting matrix with an uncorrelated propagation of the 
errors. Since the results given in these three tables are equal within errors, we chose to combine 
them to reduce their statistical errors rather than to perform a chiral extrapolation. 



1 2 3 4 5 6 7 

1 0.825(7) 0. 0. 0. 0. 0. 0. 

2 0. 0.882(38) -0.111(41) -0.009(12) 0.010(10) 0. 0. 

3 0. -0.029(69) 0.962(92) 0.013(22) -0.011(25) 0. 0. 

4 0. -0.04(12) -0.01(13) 0.924(42) -0.149(35) 0. 0. 

5 0. 0.17(18) 0.08(23) -0.042(55) 0.649(63) 0. 0. 

6 0. 0. 0. 0. 0. 0.943(8) -0.154(9) 



7 0. 0. 0. 0. 0. -0.0636(53) 0.680(11) 



The choice of external momenta specified by Eq. IA1I is non-exceptional since no partial 
sum of these momenta vanish (if their signs are chosen so that all four momenta are incoming) 
and is the choice used in Refs. |26| and [l6[. Such a choice of kinematics is expected to result 
in normalization conditions which are less sensitive to non-zero quark masses and QCD 
vacuum chiral s ym metry breaking than would be the case if an exceptional set of momenta 



had been used 
Ref. 



32|. The resulting matrix Z lat ^ RI (/i, a)/Z 2 obtained for /i = 2.15 GeV in 



is given in Tab. IVHI 

Since these RI/MOM renormalization conditions are being imposed for off-shell, gauge- 
fixed external quark lines, we must in principle include a larger number of operators than the 
minimal set of seven independent operators which can represent all gauge invariant matrix 
elements between physical states of Hy/. Therefore, we must also employ a correspondingly 
larger set of conditions to distinguish among this larger set of operators. This larger set 
of operators is required if we are to reproduce with these RI operators all the gauge-fixed, 
off-shell Green's functions that can be constructed using the original, chiral basis of lattice 
operators Q[. Thus, as stated in Sec. [V] the relations given in Eq. [22] between the seven 
lattice and the seven RI operators are valid only when those operators appear in physical 
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matrix elements between on-shell states. For this equation to be valid when the operators 
appear in the off-shell, gauge- fixed Green's that define the RI scheme, additional RI/MOM- 
normalized operators must be added. 

However, our ultimate goal is to evaluate on-shell, physical matrix elements of these 
operators. For such matrix elements there are only seven independent operators and we 
can collapse the expanded set of operators referred to above back to the seven, four-quark, 
chiral basis operators Q Rl . This is the meaning of the 7x7 matrix Z lat ^ m matrix given 
in Tab. IVII1 gauge symmetry and the equations of motion must be imposed to reduce to 
seven the Rl-normalized operators to which the seven lattice operators are equated. In the 
calculation of Z lat ^ RI presented in Ref. j^fj] such extra operators are neglected. For all but 
one, this might be justified because these operators enter only at two loops or beyond and 
the perturbative coefficients that we are using in later steps are computed at only one loop. 
A single operator, given in Eq. 146 of Ref. jg] and Eq. 12 of Ref. [l(| does appear at one loop 
but has also been neglected because it is expected to give a smaller contribution than other 
two-quark operators with quadratically divergent coefficients whose effects are indeed small. 
A final imperfection in the results presented in Tab. IVHI is that the subtraction of a third 
dimension-four, two-quark operator which contains a total derivative was not performed. 
However, the effect of subtracting this third operator is expected to be similar to those 
of the two operators which were subtracted, effects which were not visible outside of the 



statistical errors (see e.g. Tabs. XIV and XVIII in Ref. |fj]). 

In the second step we convert the seven RI operators obtained above into the MS scheme: 

Ql m =J2(l + Ar RI ^) _ Qf. (A2) 

3 13 

Here the indices % and j run over the set {1, 2, 3, 5, 6, 7, 8} corresponding to the chiral basis of 
the operators Q, defined above and a set of operators Q« with .dentinal ehiral property 

16l | . We use the computational framework described in Ref. [16] 



■i 3 

which are defined in Ref. 



and the resulting 7x7 matrix Ar RI ~^ MS is given in Tab. VIII of that reference. As in the case 
of Eq. [22J the two sets of seven RI and MS operators are related by this 7x7 matrix only 
when appearing in physical matrix elements. Since the values in this table were obtained for 
the case that the wave function renormalization constant for the quark field is the quantity 
Zg it is that factor which we use to extract Z Ut ^ Rl from the matrix Z lat ~ >RI /i?g given in 
Tab. IVHl For our f3 = 2.13, Iwasaki gauge ensembles z\ = 0.8016(3). (Note, z\ is the same 
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as the quantity Z' introduced in earlier, exceptional momentum schemes 33].) 

A third and final step is needed before we can combine the Wilson coefficients with the 
matrix elements determined in our calculation to obtain the physical amplitudes A Q and 
A 2 . The 7x7 matrix given in Tab. Vlll of Ref. gives us MS operators _defined in the 
chiral basis. However, the Wilson coefficients which are available in Ref. [l7| are defined 
for the ten operator basis referred to as basis I in Ref. |16 ]. The conversion between the 
linearly independent, seven operator basis and the conventional set of ten linearly dependent 
operators is correctly given by the application of simple Fierz identities for the case of the 
lattice and RI/MOM operators. As is explained, for example, in Ref. 16j, this procedure 
is more complex for operators defined using MS normalization. Here subtleties of defining 
7 5 in dimensions different from four, result in ten MS-normalized operators, Qf s , which are 
not related by the usual Fierz identities, with Fierz violating terms appearing at order a s . 

Thus, the conventional ten MS-normalized operators Qf 18 which appear in Eq. H]must be 
constructed, again through one-loop perturbation theory, from the seven operators Q'^: 



Q 



MS 



AT, MS ) Q' MS , 

/ in 



(A3) 



in the notation of Ref. 



The 10 x 7 matrices, T and AT 7 MS are given in Eqs. 59 and 65 of 
that reference. (The subscript / on the matrix AT^ S identifies the particular ten-operator, 
MS basis required by the Wilson coefficients of Ref. 

This entire set of non-perturbative and perturbative transformations can be summarized 
by the following equation which expresses the ten MS-normalized operators Qf 1 in terms 
of the seven, chiral basis, lattice operators whose matrix elements we actually compute: 



Q 



MS 



E 

j 

E 



t + at; 



10x7 



1 + Ar R I_> s ) C^iat->Ri\ 

/ 7x7 



7x7 



Qj 



lat 

3 



2-lat->Ms\ 

/ 10x7 



(A4) 
(A5) 



'j 



where the subscripts indicate the dimensions of the matrices being multiplied and the matrix 
Z)f^ m is used in Eq. [25] 

The physical matrix elements listed in Tabs. [Ill and [V] are obtained by using Eq. IA5I to 
determine the matrix elements of the ten conventional operators in term of the matrix 
elements of the seven lattice operators Qj. These ten matrix elements are then combined 
with the twenty Wilson coefficients computed for the renormalization scale n — 2.15 GeV 
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TABLE VIII: Wilson Coefficients in the MS scheme, at energy scale [i = 2.15GeV. 



i 






1 





-0.29829 


2 





1.14439 


3 


0.024141 


-0.00243827 


4 


-0.058121 


0.00995157 


5 


0.0102484 


-0.00110544 


6 


-0.069971 


0.00657457 


7 


-0.000211182 


0.0000701587 


8 


0.000779244 


-0.0000899541 


9 


-0.0106787 


0.0000150176 


10 


0.0029815 


0.0000656482 



using the formulae in Ref. [13]. The values obtained for these Wilson coefficients are listed 
in Tab. IVHT1 

Note, there are many important details of the RI/MOM renormalization procedure, such 
as the subtraction of dimension three and four operators, which are not repeated here because 
they are already discussed with some care in Refs. 6] and [if]]. 
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